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Introduction 


In  the  past  few  years  the  description  of  the  earth’s  gravity  potential 
in  terms  of  spherical  harmonic  coefficients  has  been  extended  to  degree 
180  and  in  special  cases  to  higher  degrees  (Rapp,  1978,  1981),  and  Lerch 
et  al.  ( 1981) .  These  high  degree  expansions  can  be  used  to  evaluate  quantities 
such  as  geoid  undulations,  height  anomalies,  gravity  anomalies,  gravity 
disturbances,  deflections  of  the  vertical,  etc.  To  do  this  efficient  com¬ 
puter  programs  are  needed.  The  purpose  of  this  report  is  to  describe  one 
Fortran  computer  program  that  can  be  used  for  these  calculations. 


Theory— Basic  Equations 

The  gravitational  potential,  V  ,  in  spherical  harmonics  can  be  written 
as: 

<f>"  j. 

n=£  m=0 

+  lm  sinmX)  Pnm  (sin*)]  (1) 

where:  kM  is  the  geocentric  gravitational  constant; 

r  is  the  geocentric  radius; 

is  the  geocentric  latitude; 

X  is  the  "geocentric"  longitude; 

Cnm»5nm  are  the  fully  normalized  potential  coefficients; 

a  is  the  scaling  factor  associated  with  the  coefficients. 

The  disturbing  potential,  T  ,  is  the  difference  between  the  actual  potential 
(V)  at  a  point  and  the  "normal"  potential  at  the  corresponding  point.  For 
our  purpose  the  normal  potential  will  be  that  associated  with  an  equipoten- 
tial  reference  ellipsoid  of  defined  parameters.  We  have: 

T(r ,^,X)  =  V(r,ip,X)  -  U(r ,ip,A)  (2) 


-1- 


The  potential  associated  with  U  can  be  described  by  an  even  degree  zonal 
harmonic  expansion.  We  can  write: 


T(r,*.x> .  (sa^ad  ♦  a.  i  (£)»  I 

r  r  n-2  r  m-0 

(C^cosmx  ♦  S,,.  stn.x)  (3) 

where  kMg  is  the  mass  of  the  reference  ellipsoid  and  C*m  are  the  dif¬ 
ferences  between  the  actual  coefficients  and  those  implied  by  the  reference 
equipotential  ellipsoid.  We  have: 

«?.,■ «...  - 

c...  -  c...(ref)  (4> 

C*  *  C  -  C  (ref) 

6  »0  6  >D  6  »0  '  ' 

In  most  cases  we  assume  kM  is  equal  to  kM^  so  that  (3)  becomes: 

CO 

T(r,ipfX)  •y-  l  (f)n  l  (C*m  cosmX  +  §nm  sinmX)  Pnm(sin^)  (5) 
n*2  m*0 

In  classical  gravimetric  geodesy  we  discuss  geoid  undulations,  N  , 
and  geop-spherop  separations,  Nr  .  If  W0  is  the  potential  of  the  geoid 
and  Uo  is  the  potential  on  the  surface  of  the  reference  ellipsoid  the 
geop-spherop  separation  is  (Heiskanen  and  Moritz,  1967,  Section  2-19): 

N(r,*,x)  -  (6) 

4 

where  y  is  normal  gravity.  In  most  cases  we  take  WQ *  Uo  so  that  (5) 
becomes: 


N(r,if»,X) 


IkjM J 

Y(r,ip) 


The  non-classical  procedure  uses  the  concept  of  the  disturbance  poten¬ 
tial  at  some  surface  points  and  introduces  the  term  height  anomaly,  q  . 


Let  W(r,i/>,A)  define  the  gravity  potential  and  U(r,tp,A)  the  normal  gravity 
potential  at  the  same  point.  Then: 

T(r,!j>,A)  =  W(r,4>,A)  -  U(r,^,A)  (8) 

We  can  introduce  the  geopotential  number,  Cp  ,  with  respect  to  a  reference 
potential,  Wo  ,  such  that 

Cp  ■  W(r,iM)  -  Wo  (9) 

The  normal  height  of  P  ,  H*  ,  can  be  computed  from  Cp  (Heiskanen  and 
Moritz,  1967,  section  4.5).  Letting  h  be  the  geometric  height  of  P 
above  the  reference  ellipsoid  the  height  anomaly  is: 

C  =  h  -  H*  (10) 

In  terms  of  the  disturbing  potential  we  can  write: 


This  equation  is  the  same  as  (7)  but  there  will  be  a  conceptual  (but  small) 
difference  when  comparing  normal  heights,  height  anomalies,  geoid  undulations 
(N)  and  orthometric  heights  H  .  Specifically  we  have  (ibid,  section  8-12): 

h»H+N»H*  +  5  (12) 

For  our  purposes  we  consider  the  disturbing  potential  to  be  given  by 
equation  (5)  with  the  calculation  of  the  height  anomaly  by  (11).  For  the 
calculation  of  geoid  undulations  we  would  also  use  equation  (11)  but  with 
the  evaluation  of  T  on  the  geoid  by  the  appropriate  choice  of  r  .  Although 
the  convergence  of  the  Infinite  series  for  T  is  a  formal  concern  the  cal¬ 
culations  of  Jekell  (1981)  with  high  degree  finite  series  show  that  there 
is  no  practical  concern. 


The  gravity  anomaly  is  a  vector  that  can  be  expressed  in  the  classical 
and  pon-classical  forms.  In  eithr"  case  the  general  relationship  is  the 
same  between  the  disturbing  potent  1  and  the  anomaly  although  there  is 
a  conceptual  difference.  For  the  anomaly  component  in  the  vertical  direction 
(h)  we  have  (Heiskanen  and  Moritz,  p.  967,  p.  84,  and  298): 

Agh(r,M)  *  -  yj-  +  7-  57-  T(r,tM)  (13) 


For  the  classical  anomaly  at  the  geoid,  T  is  evaluated  there,  while  for 
the  surface  anomaly  T  is  evaluated  at  the  surface  point.  We  can  obtain 
the  radial  component  of  the  anomaly  by  writing  (13)  in  the  form: 

Agr(r,iM)  ="37-  +  7"  37"  T(r»'M)  (14) 


With  a  spherical  approximation  we  have: 

lil  .  _2_ 

Y  dr  r 

so  that  (14)  becomes 


(15) 


Agr(r,ip,A)  =  -  T(r,ip,A)  (16) 

If  we  now  take  equation  (5)  for  T  we  have: 

CO  I'j 

Ag  (r,if),A)  *  l  (n-1)  (|*)n  {  (C  cosmA 
r2  n=2  r  m=0  m 

+  §nm  sinmA>  ~Pm  (sin*>  ™ 

The  deflection  of  the  vertical  represents  the  angular  difference  between 
the  normals  to  the  actual  gravity  equipotential  surface  and  the  normal  equi- 
potential  surface.  For  a  deflection  in  an  arbitrary  direction  (s)  we  can 
write  (Pick  et  al.  ,  1973,  p.  257): 

a  -  I  'LL 

9  gas 
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(18) 


where  s  lies  in  the  plane  tangent  to  the  normal  equipotential  surface. 
Normally  the  total  deflections  is  expressed  in  a  meridian  (?)  component 
and  a  prime  vertical  (n)  component.  In  the  meridian  we  have,  with  sufficient 
accuracy  ds  =  rd^>  ,  and  in  the  prime  vertical,  ds  =  rcos^dX  .  Thus  the 
deflections  of  the  vertical  are: 


±11  _  l  jrr 

gr  dip  ’  n  grcosy;  3X 


(19) 


As  pointed  out  in  Pick  et  als  (1973,  p.  307)  the  derivatives  and 

"are  the  derivatives  of  the  disturbing  potential  with  respect  to  the 
appropriate  direction  assuming  that  H  and  X  ,  and  H  and  >t>  ,  respectively, 
are  constant".  (H  corresponds  to  height  and  <p  latitude.)  Thus  it  is 
possible  to  use  (5)  for  T  to  obtain  the  deflections.  We  then  have,  letting 
9p  =  Y  (r.ijj): 


5  *  -  fr  l  <f>"  <C;m  d 

i '  n -i 


(20) 


n  * 


kM  l  (|-)n  l  m(C*  (-sinmX)  +S  cosmX)  Pnm(s1ni|>)  (21) 


Yr  cosi^n=2 


m*0 


To  obtain  the  deflections  in  seconds  multiply  the  above  equations  by  the 
radian  conversion  factor. 


The  gravity  disturbance  vector  is  defined  as  the  difference  between 
gravity  at  a  point  and  normal  gravity  at  the  same  point.  We  have  (Heis- 
kanen  and  Moritz,  1967,  p.  84): 


3  =  Op  -  YP  =  grad  T 


(22) 


The  radial  component  of  the  gravity  disturbance  can  be  defined  as  (ibid, 
p.  85) 


,  -  9T 
r  3r 


(23) 


In  some  cases  (ibid,  p.  233)  the  minus  sign  is  not  used  but  we  retain  (23) 
as  our  defining  equation.  Noting  that  (23)  appears  in  (16)  we  can  avoid 


a  direct  evaluation  of  6r  by  computing  it  from  (16)  after  Ag  and  I 
(or  N ( c) )  have  been  computed.  We  have: 


Sr  =  Ag  +  |-  T 
r  r 


(24) 


The  other  two  components  of  the  gravity  disturbance  vector  are  defined  as 
(ibid,  p.  285) 


p  _  _1_  3T  x  =  1  ST 

°ip  ~  r  dip  ’  °A  r  cos^  3X 

Comparing  these  quantities  to  (19)  we  see 


(25) 


(26) 


where  we  have  let  g  =  y  .  Thus  once  £  and  n  are  computed  it  is  a  simple 
matter  to  calculate  the  two  disturbance  components  6^  ,  6^  . 


In  summary  we  are  given  a  set  of  fully  normalized  potential  coefficients. 
From  these  coefficients  we  remove  the  values  implied  by  an  equipotential 
reference  ellipsoid.  This  leaves  us  with  the  expression  for  the  disturbing 
potential  T  (equation  5).  We  then  can  compute  height  anomalies  from  equation 
(11).  gravity  anomalies  from  equation  (17),  the  deflections  of  the  vertical 
from  (20)  and  (21),  and  the  radial  gravity  disturbance  from  equation  (24). 


Theory— Auxilary  Relationships 

To  implement  the  equations  discussed  in  the  previous  section  a  number 
of  additional  quantities  are  needed.  These  are  now  discussed. 

The  Reference  Potential  Coefficients 

Given  four  parameters  defining  an  equipotential  reference  ellipsoid 
all  the  even  degree  zonal  harmonics  are  explicitly  defined.  For  our  program 
it  is  sufficient  to  use  only  the  zonal  terms  to  degree  six  as  taken  from 
Cook  (1959).  We  have  given: 


-6- 


a  *  equatorial  radius 
kM  =  geocentric  gravitational  constant 
u  =  angular  velocity 
f  =  flattening. 


Then  compute  m  : 


The  zonal  coefficients  in  the  Jn  form  are: 

J2  3  y[f(l-if)  -  im(l-  yf  +  f2)] 

Ju  =  -  ^fd-if)  (7f(l-if)  -  5m(l  -  yf ) ) 
Je  =  f'  (6f  -  5m) 


(27) 

(28) 

(29) 

(30) 


These  coefficients  are  related  to  the  fully  normalized  C  coefficients 
through  the  following: 


'n0 


/2n+T 


(31) 


Calculation  of  ,  and  r 

Generally  the  latitude  point  will  be  specified  as  a  geodetic  latitude. 
Formally  this  latitude  should  be  with  respect  to  an  ellipsoid  whose  center 
is  at  the  center  of  mass  of  the  earth.  The  geodetic  latitude  must  be  converted 
to  a  geocentric  latitude  and  the  geocentric  radius  must  be  computed.  Given 
0,  X,  and  h  the  rectangular  coordinates  of  the  point  are  (Rapp,  1981  equation 
60). 


X  *  (N+h)  cos<j>  cosX 
Y  =  (N+h)  cos<J>  sinx 
Z  *  (N ( 1  -e  2)  +  h)  sin<p 
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(32) 


where  N  is  the  prime  vertical  radius  of  curvature: 


(1  -  e2  sin2^)^ 

The  geocentric  radius  is  then: 

r  =  (X2  +  V2  +  Z2)i 

The  geocentric  latitude  is  then 

V  =  tan'1  1  .  (35) 

/X2  +  Y2 

Calculation  of  y 


(33) 


(34) 


Normal  gravity  is  needed  in  the  evaluation  of  the  height  anomaly  and 
deflections  of  the  vertical.  A  high  degree  of  accuracy  is  not  needed  *or 
this  calculation  as  the  number  of  digits  in  the  final  quantities  is  usually 
only  two  to  four.  In  our  case  we  choose  to  evaluate  normal  gravity  for 
the  point  on  the  ellipsoid  and  then  modify  this  value  in  a  linear  fashion 
for  the  height  of  the  point  above  the  ellipsoid.  The  normal  gravity  on 
the  ellipsoid  is: 


I  +  k  sirr2^ 
b  A  -  e 2  s  i  n 2  <i> 

The  value  of  y  at  the  height  h  above  the  ellipsoid  is: 


(36) 


Yh  =  y  -  0.3086  x  10"5  h 
where  y  is  in  meters/s2  and  h  is  in  meters. 


(37) 


Calculation  of  Pnm  and  Its  Derivative 

The  generation  of  the  fully  normalized  associated  Legendre  functions 
and  its  first  derivative  is  critical  to  any  calculation  involving  spherical 
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harmonic  expansions.  In  choosing  an  algorithm  one  must  consider  the  speed 
and  the  stability  and  accuracy  of  the  procedure.  In  the  past  few  years  a 
number  of  different  equation  sets  have  been  described  in  the  literature. 

For  this  program  we  have  chosen  subroutine  LE6FDN  that  is  described 
by  Colombo  (1981,  p.  131).  Colombo  has  carried  out  a  number  of  tests  to 
investigate  the  stability  of  the  equations. 

The  subroutine  is  written  such  that  the  needed  functions  for  a  given 
order  m  and  all  degrees  to  the  highest  maximum  degree  are  computed  in 
one  call  to  the  subroutine.  The  subroutine  is  repeatedly  called  for 
0  s  m  s  N  where  N  is  the  maximum  degree  being  used  in  the  expansion. 

For  discussion  purposes  visualize  the  associated  Legendre  functions 
in  a  lower  triangular  matrix  where  the  rows  correspond  to  degree  n  and 
the  columns  correspond  to  order  m  . 

For  a  given  m  ,  the  subroutine  first  calculates  for  0  <  r.  £  m  the 
diagonal  elements  corresponding  to  the  diagonal  passing  through  the  n  =  m 
location.  We  have: 

Pnn.(cose>  *  sfne  pn-l,n-lCcos0) 

Poo(cose)  =  1.0  (38) 

Pii(cose)  =  sine 

Then  the  following  element  is  computed: 

Pn+i#n  (cose)  =  /2n+3  cose  Pn>n  (cose)  (39) 

with  n  =  m  .  Then  the  following  recursive  relationship  is  used  to  calculate 
the  remaining  values  of  P  for  m+2  s  n  t:  N  . 


Pnm<c°s®>  *  /in-m)i^)‘)  c°se  •Vi,.  <c°se’ 


n  s  2  ,  (n-2)  >  m  *  O' 


Note  that  3  is  the  polar  angle  given  by 


(40) 


3  =  90°  -  Ip  (41) 

Singh  (1982)  has  pointed  out  ways  to  improve  the  calculation  of  the 
associated  Legnedre  functions  by  applying  a  scaling  factor,  such  as  1  x  1072 
in  the  recursive  procedure.  This  procedure  was  tested  in  actual  calcula¬ 
tions  of  £  ,  Ag  etc.  No  difference  in  numberical  values  was  seen  when 
the  scale  factor  was  and  was  not  applied.  Consequently  we  did  not  implement 
the  scaling  operation.  Doing  so  might  avoid  some  underflow  messages,  but 
would  not  changes  results  when  using  expansions  to  degree  180. 


For  the  £  component  calculations  we  need  the  derivative  of  P  .  Colombo 
(1981)  implemented  the  following  procedure: 


dPnn  (cose) 

de 


-  r  2n+lii 

’  L“2fT] 


(sine 


dP 


n-l.flri 

de 


+  cose 


P 


n-l,n-l 


(cose)) 


(42) 


After  these  values  are  computed  for  a  given  m  up  to  a  given  N  ,  then 
we  have: 


dP 


=  (sine)_1(n  Pnjn  (cose)  cose 
-  F*  Pn- 1  ,m  <«*•» 


The  starting  value  is 


dP  o  o 

de 


=  0 


(43) 


(44) 


Due  to  the  occurrence  of  (sine)-1  this  subroutine  can  not  calculate  the 
derivatives  at  the  poles. 

Since  we  want  the  derivative  of  P  with  respect  to  ip  we  note  that: 


dP  dP 

dijj  ~  "  "3T 


(45) 


The  calculation  of  sinmA  and  cosm^ 

The  generation  of  sinmA  and  cosmA  is  done  through  the  following 
recursion  relationships: 

sinmA  =  2cos  A  sin(m-l)  A  -sin(m-2)  A 

cosmA  =  2cos  A  cos(m-l)  A -cos(m-2)  A  (46) 

These  relationships  are  useful  for  point  calculations  but  are  inefficient 
for  use  if  a  set  of  points  at  a  uniform  longitude  interval  are  being  used. 


Geodetic  Constants 


For  the  evaluation  of  the  reference  potential  coefficients,  the  geocen¬ 
tric  radius  vector  etc,  we  need  to  adopt  a  set  of  constants.  We  used  the 
values  of  the  Geodetic  Reference  System  1980.  We  have: 

a  =  6378137  meters 

kM  =  3986005  x  10 8  m3  s"2 
oo  =  7292115  x  10"11  rad  s"1 
e2  =  0.006  694  380  022  90 

f  =  0.003  352  810  681  18 

y  =9.780  326  7715  ms'2 

k'  =  0.001  931  851  353 

These  constants  are  used  in  the  calculation  of  the  reference  potential  coef¬ 
ficients  (for  a  flattening  that  is  read  into  the  program),  the  geocentric 
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radius,  and  normal  gravity.  These  constants  can  easily  be  changed  in  the 
program. 

It  is  critical  to  note  that  the  use  of  the  above  constants  does  not 
mean  that  the  geoid  undulation  (for  example)  refers  to  the  GRS80  reference 
ellipsoid.  This  is  because  the  zero  order  term  in  T  has  been  set  to  zero. 
The  real  reference  ellipsoid  is  that  one  which  best  fits  the  geoid  and  this 
may  or  may  not  be  GRS80. 


The  Program 

The  program  written  to  implement  the  equations  previously  described 
is  given  in  Appendix  A.  This  Fortran  program  was  run  on  an  Amdahl  470  V-8 
machine  using  double  precision  computations. 

The  program  is  currently  designed  for  point  by  point  calculation.  In 
this  case  the  input  information  is  as  follows: 

1.  NMAX,  F  (I3.F10.4) 

NMAX  is  the  highest  degree  to  be  used  in  the  expansion,  F  is  1/f  which 
is  the  inverse  flattening  of  the  reference  ellipsoid  to  which  the  computed 
quantities  are  to  be  referred. 

2.  The  fully  normalized  potential  coefficients  are  read  from  tape  or  disk 
file  in  the  form  of  (n,m,Cnm,5nm) .  The  arrangement  of  the  input  is  in  order 
of  degree,  i.e.  from  lowest  to  highest  degree.  However  the  storage  location 
for  the  coefficients  is  computed  from  the  given  n  and  m  values.  In 
this  program  all  coefficiecns  are  stored  in  double  precision.  Space  can 

be  sa'  u  iy  storing  in  single  precision. 

3.  The  coordinates  of  the  points  at  which  Q  and  the  other  quantities  are 
to  be  computed.  Specifically  ($,A,  h)  where  $  is  the  geodetic  latitude, 

A  is  the  longitude  and  h  is  the  height  above  the  reference  ellipsoid. 

The  current  format  is  (3F10.1).  The  last  point  is  signaled  by  an  end  of 
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file  (/*)  card.  Points  having  the  same  latitude  should  be  grouped  together 
as  in  this  case  the  associated  Legendre  functions  and  their  derivatives 
are  not  re-computed. 

The  output  is  printed  across  the  page  under  column  headings:  LAT, 

LON,  HEIGHT,  UNDU,  ANOM,  DIST,  XI,  ETA.  Although  the  output  is  given  two 
decimal  digits,  actual  accuracy  is  considerably  poorer  than  this  because 
of  the  errors  in  the  potential  coefficients. 

The  values  computed  by  this  program  have  been  checked  against  another 
program  written  by  Tscherning  and  Goad  (1982,  private  communication).  All 
values  checked  agreed  to  two  decimal  digits. 

The  computer  time  needed  for  a  single  point  calculation  (after  the 

potential  coefficients  are  input)  is  0.46  seconds  with  an  expansion  complete 
to  degree  180  on  an  Amdahl  470  V/8.  A  calculation  of  points  on  a  12°x  12° 
grid  at  1°  intersection  took  21.9  seconds.  If  a  limited  grid  of  undulations 
or  anomalies  are  to  be  generated  the  program  described  by  Rizos  is  the  most 
efficient  procedure  to  date.  If  a  global  grid  is  being  generated  the  program 
SSYNTH  described  in  Colombo  (1981)  is  the  most  efficient. 

For  checking  the  results  of  the  program  the  values  of  C  ,  Ag  ,  5  , 

5  ,  n  have  been  computed  at  five  test  points  using  three  different  sets 
of  potential  coefficients  to  degree  180.  These  values  have  been  computed 
with  respect  to  an  ellipsoid  which  has  the  flattening  of  GRS80  and  are 
given  in  Table  1. 
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Table  1 


Sample  Computed  Values 
(reference  flattening  =  1/298.257222) 


<t>° 

X" 

h(m) 

C  (m) 

Ag(mgals)  {(mgals) 

c 

n" 

21° 

1° 

0 

Rapp78 

34.46 

20.07 

30.65 

0.68 

-0.25 

Rapp81 

30.56 

7.73 

17.11 

0.60 

0.40 

GEM IOC 

28.37 

4.12 

12.83 

-0.10 

0.21 

2  r 

45° 

0 

Rapp78 

-11.19 

-4.75 

-8.19 

-5.41 

11.12 

Rapp81 

-9.58 

-5.55 

-8.49 

-4.24 

10.63 

GEM10C 

-9.68 

-8.46 

-11.43 

-2.23 

8.98 

5° 

79° 

0 

Rapp78 

-104.42 

-84.60 

-116.62 

-1.43 

0.64 

Rapp81 

-107.48 

-91.84 

-124.81 

0.02 

0.65 

GEM10C 

-106.20 

-87.66 

-120.23 

-1.13 

-1.04 

5° 

79° 

10000 

Rapp78 

-103.58 

-78.90 

-110.52 

-1.63 

0.35 

Rapp81 

-106.58 

-85.49 

-118.02 

-0.22 

0.50 

GEM10C 

-105.35 

-80.51 

-112.66 

-1.12 

-0.93 

00 

"*4 

0 

21° 

0 

Rapp78 

15.43 

-1.46 

3.32 

1.32 

2.37 

Rapp81 

20.23 

8.86 

15.12 

0.81 

1.86 

GEM10C 

18.38 

3.58 

9.26 

2.59 

4.05 

Summary 


This  report  describes  a  Fortran  computer  program  that  can  be  used 
for  the  calculation  of  c  ,  Ag  ,  6  ,  £  ,  n  which  are  dependent  on 

a  set  of  fully  normalized  potential  coefficients.  The  program  has  been 
set  to  work  to  degree  180  and  it  can  be  extended  higher. 

The  equations  used  for  the  calculations  are  to  some  extent  spherical 
approximations.  However  literal  interpretation  of  certain  quantities  would 
be  formally  correct  (e.g.  the  radial  component  of  the  gravity  disturbance). 
Correction  terms  for  spherical  harmonic  expansions  evaluated  considering 
the  ellipticity  of  the  earth  are  described  to  some  extent,  in  Jekeli  (1981, 
section  4). 

The  input  quantities  to  the  program  are  geodetic  latitude,  longitude 
and  height  above  the  ellipsoid.  In  theory  these  quantities  should  be  given 
with  respect  to  a  geocentric  ellipsoid.  In  practice  the  use  of  non-geocentric 
coordinates  would  cause  small  but  systematic  errors  in  the  results. 

The  computed  quantities  refer  to  a  geocentric  ellipsoid  whose  flattening 
is  an  input  parameter.  The  size  of  this  ellipsoid  is  not  specifically 
defined  because  the  zero  degree  term  in  the  disturbing  potential  expansion 
has  been  set  to  zero.  In  most  applications  the  equatorial  radius  of  the 
ellipsoid  is  the  current  best  estimate. 

The  computer  program  of  this  report  has  been  checked  against  other 
programs  with  excellent  agreement.  The  stability  of  the  algorithms  for 
the  associated  Legendre  functions  has  been  checked  by  Colombo  (1981)  and 
by  Singh  (1982).  For  some  applications  at  high  latitude  underflows  may 
occur  in  the  computations.  These  are  machine  dependent  quantities  and 
can  be  turned  off  if  desired. 

Other  procedures  have  been  developed  that  extend  the  derivatives 
of  the  potential  to  the  second  derivative  (Tscherning  and  Poder,  1981). 

In  addition,  problems  at  the  pole  that  exist  with  our  current  program  (for 
the  derivative  of  Pnm)  are  avoided  with  the  Tscherning/Poder  application 
of  the  Clenshaw  summation. 
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//  JCB  • XXXXX.XXXXXXXXX' ,'RAPP,R.H.'  , 

//  TIMt=( 0,40) |REGIGN«1024K 
/* J08PARM  LINES=J*0C,DISKI0=2400,V=R 
//SI  EXEC  FORTQCG 
//FORT .SYS  IN  DO  * 

C  THIS  PROGRAM  WAS  PUT  IN  ITS  PRESENT  FORM  BY  R.H.  RAPP  IN  AUG  L9b2 
C  THE  PROGRAM  IS  A  MODIFICATION  CF  PROGRAM  F379 
C  PUINT  COMPUTATION  FROM  HARMONIC  COEFFICIENTS 

C  DIMENSIONS  CF  P,Q,hC,HS  MUST  BE  AT  LEAST  I  {  MA  XN+ 1)  *(  •«A  XN*  2  ))  U  , 

C  DIMENSIONS  OF  SINML ,CCSML , SCRAP  MUST  Be  AT  LEAST  MAXN, 

C  WHERE  MAXN  IS  MAXIMUM  ORDER  OF  CuMPuTATluN 

C  THE  CURRENT  DIMENSIONS  ARE  SET  FUR  A  MAXIMUM  DEGREE  OF  loJ 
IMPLICIT  REAL*8  (A-H,G-Z) 

REAL *8  P l 16471) , SCRAP t 131 ) ,RLEG( 131 ) , DLEG( 131 ) , RLNNt 181) , 
*P0ER(la4  7l),SINMLU81)  ,CCSML(  181  ) 

REAL*8  HC(  16471)  ,HS(  16471) 

DATA  RAD/57. 2957795130823200/ 

C  F  IS  THfc  REFERENCE  INVERSE  FLA  I  TEN  I NG» NM  AX  IS  THE  MAXIMO*  DEGREE 
010  READ ( 5 ,  900)  NMAX  »F 
900  FORMAT ( I  3 ,F10.4 ) 

WRITEl6,910)  NMAX , F 

910  FGRMAT l 1H1 »// •  MAXIMUM  DEGREE  = • , 1 4 , BOX , • 1 / F  =  • , F 6 . 4/// ) 
F=1.0D0/F 

CALL  OHCS INI NMAX, F ,RJ2  »R J4,RJ6,HC ,HS) 

C  SETTING  IFLAG-=0  FORCES  LEGEnORE  FJNCTION  DbnIVATIVtS  TC  3E  TAKEN 
IFLAG=0 
IR  =0 

K=NMAX+1 

FLATL=90.0U0 

978  FO^MA^f IhI?/ /  '  LAT  »,*  LCN  HEIGHT  JNu., 

*•  ANOM  OIST  XI  ETA') 

C  READ  GEODETIC  LA  T  I  T'JDE  ,  LONGI  T'JDE  ,  AND  HEIGHT  ABOVE  THE  ELLIPSOID 
030  RE AO (5,930, £  ND*G90 )  FL AT , FLON , H T 
930  F0RMATC3F10.il 

C  COMPUTE  THE  GEC'CENTRIC  LAT  I T  U9E »  GCOCENTk  I C  RADIUS  ,  NCRMmL  GRAVITY 
CALL  RAOGRAtELAT ,FLON,HT ,RLAT,GR,RE I 
IFlFLATL.Ew.rLAT)  GO  Tu  040 
R  L  AT  1 =RLA  T 

RLAT»l.57u  79  6326  79 4 8 96 600- RL AT 

FLATL*FLAT 

DO  25  J«l,K 

M«J-i 

CALL  LtGFDNl M ,RL  AT ,R LEG, OLEG ,NMAX , I R , KLNN , I  FLAG ) 

DO  26  I  * J  t  K 
N=  I- 1 

L0C=(N*(N+1) )/2>M+i 
POER (LOCI ®DLEG( I ) 

26  P( LOC I aRLEGC I ) 

25  CONTINUE 
040  kLUN*FLCN/RAO 

CALLDSCML  ( R LGN , NM AX  »  S I NML , C03ML ) 

CALL  HUNDU  ( U  » DC , 0 1 S  »  X  I  , E  T A  , HT  ,  NMAX ,  P  ,  PEER,  hC  ,  HS  ,  G  I NML  ,  CgML  , 
♦GR.RE.RLAT1) 

WRITE! 6,940)  FL AT*FLON , HT *U*DG *DI S. XI , ET  A 
940  FORMAT ( 2E9.4,IF09.2,5F10.2) 

GO  TO  30 
90  STUP 
END 

SUBROUTINE  HtJNDU  IUNDU ,  ANOM,  D  I  ST,  X  I  ,  ET  A  ,HT  ,NMAX  ,  P ,  PULR  ,HC  ,HS, 

♦  SINML  »CGSML ♦ GR , RE , ANG) 

IMPLICIT  REA L*3  ( A-H ,Q-Z ) 

DIMENSION  SINML < 1) ,CUSML( I) ,P( l) ,POER ( l) 

REAL *8  HC(1)  ,HSU) 
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C  CONSTANTS  FOR  GRS80 

DATA  GM/.  39B  600  5D1  3/,  AE/6  376  137.000/,?.  HO/206264. 80600/ 

AKxAE/RE 

ARN=AK 

A=Q  .  0 

3=0.0 

XI=0.0 

ETA=C.O 

K=3 

DO  030  N* 2, NMAX 
ARN=ARN*AR 
K  =  K*  1 

SUM  =  P ( K ) *HC i K ) 

SUMl  =  PQ£RU)*riC{  K) 

SUM2=0. C 
00  020  M=  L ,  N 
K=K  + 1 

TFMP=HC(K)*CCSML  ( .*1 )  *HS IK)  * S I NML CM) 

SUMl=SuMl+PDER(K)*TEMP 

SUME=SUM2  +P I K  J  *M*I -HC  I  K  i  *S I  NML l M  J +HS < K  J *COSML  i  M ) ) 

020  SUM=SUM*  P  { K )  *TEMP 

B=B*  SUM* ARN* ( N- 1 i 
X  I  *  X  l  ♦  S  l»M  1  *  A  kN 
ETA=ETA*SUM2*ARN 
30  A=A*S'JM*ARN 

UNOU=A*GM/(GR*RE J 
ANCM=8*GM/RC**2*L.D3 
0IST=ANCM+2. *UNDU*U«/R£*1.05 

C  THE  SIGN  OF  *  IN  XI  OCCURS  DUE  TO  THE  DERIVATIVE  3EING 
C  WITH  RESPECT  TO  POLAR  DISTANCE  NCT  LATITUDE 
XI=*RHO*GM/ (G«<*RE**2 ) *XI 
ETA=-RHO*GM/ (GR*RE**2*DC0S(ANG) )*6TA 
C  THE  UNITS  OF  THE  UNJULATIUN  ARE  METERS 
C  THE  UMTS  OF  THE  ANOMALY  ANO  DISTURBANCE  ARE  MGALS 
C  THE  UNITS  UF  THt  DEFLECTIONS  ARE  SECONOS 
RETURN 
END 

SUBROUTINE  OSCML  < RLUN . NMAX , S I NML , CCSML ) 

IMPLICIT  KEA L *d  (A-H,0-Zi 
DIMENSION  S1NMU  I)  .CUSMUl) 

A*OS INI RLUN) 

B=UCCS ( RLON ) 

S I NML  C 1 )  —  A 
CCSML (  I  )  =8 
S I NML ( 2 ) =2 • 3*6*A 
C0SML(2)=2.U*  )*B-l.O 
DO  010  M* 3, NMAX 

SINML(M)«2.0*B*SINML( M— l ) -S I NML I M— 2 l 
010  COS'«UM)=2.G*'2*COSML{M-l  I -COSML ( M-2 J 
RETURN 
ENO 

SUBROUTINE  OHCSIN  I NMAX*  F  »  J2  ,  J4 , J6 , HC ,  HS  J 
IMPLICIT  RtAL*8  (A-H,0-Z) 

R E AL  *8  J2 f  J4  *  j6 
REAL*8  HC(l» ,HS(LI 

C  THIS  VERSION  USES  IMPROVED  J2,J4,Jb  F.iOM  uCU*  PAPER!  L  E-  5  3  ) 

C  CONSTANTS  FROM  GPS  30 

DATA  FKM, CM, A/ 3. S86005D1 4»7.292 I l 50-5 , 63 78 l 3 7. ODD/ 
FM»uM**2* A**  3* ( 1 • 000— F ) /FKM 
M=  UNMAX*  I  J*(NMAX*2)  )/2 
00  DO  1  N*i,M 
HC I  NT  =0.0 
001  HSIN)*0.0 

02  READ! 12,EN0=3)N,M,C,S 

I F I N.G T .NMAX )  GO  TO  0U3 
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N* 1 N* l  N  +  l }  J/2+M  +  1 
HCtNUC 
HS(N)=S 
GU  TO  002 

3  J2=2.0D0/3.CD0*(  F*{ 1 .GDO-F/2 .000 ) -FM/ 2 . ODO* ( 1 .0Du-2 . 000/ 7 . JD0*F 

*+11.0D0*F*F/49.0D0)  ) 

J4=-4. 000/35. 0U0*F* I  1 . GD0-F/ 2 . 000 ) * l 7.0D0*F*( l . JjO-F/2 . 000 ) 
*-5.0D0*FM*l i. 000-2. 0D0*F/7.0D0J  ) 

J6=4.*F**2*(  to.*F-5.*F,M) /21 . 

HC(4)=HC(4)+J2/USCRT(5. DO)  I 

HC ( 11 )=HCt  L  L ) +04/3. 000 

HC( 22)=hC(22)+J6/0SgRT( 13.00) 

RETURN 

EN0 

SUBROUTINE  LBGFDNIM, THET A, R LEG ,GLEG , NMX *  IK »  RlNN  * [FLAG ) 

THIS  SUBROUTINE  COMPUTES  ALL  NORMALIZED  LEGENDRE  FUNCTIONS 
IN  "RLEG"  AND  THEIR  DERIVATIVES  IN  "ULEG".  ORDER  IS  ALWAYS 
M  ,  AND  COL  AT  I  TUl)E  IS  ALWAYS  THETA  (RADIANS).  VAX  I  MUM  DEGPF 
IS  NMX  .  ALL  CALCULATIONS  IN  DOUBLE  PRECISION. 

IR  MUST  BE  SET  TO  ZERG  BEFORE  THE  FUST  CALL  )C  THIS  Sud. 
THE  DIMENSIONS  OF  ARRAYS  RLEG »  JLEG »  AND  RLNN  MOST  BE 
AT  LEAST  EQUAL  TO  NMX+1  . 

THIS  PROGRAM  DOES  NOT  COMPUTE  DERIVATIVES  AT  THE  PULES  . 

IF  I  FLAG  =  L  ,  ONLY  THE  LEGtNORE  FUNCTIONS  ARE 
COMPUTED. 

ORIGINAL  PROGRAMMER  :OSCAR  L.  COLOMBO ,  DEPT.  OF  GEODLTIC  SCIENCE, 
THE  OHIO  STATE  UNIVERSITY,  AUGUST  19C0  .  ***** ********* *v ** 


IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  RLEG (II ,DL£G( 1 ) , RLNN ( 1 ) 

2,  DRTS (  1300)  ,DIRT( 13C0  ) 

NMX I  =  NMX+1 
NMX2P  =  2*NMX+ 1 
Ml  =  M+l 
M2  =  M+2 
M3  =  M+3 

I F ( IR.EQ.l)  GO  TO  10 
IR  *  1 

DU  5  N  ~  l , NMX 2P 

DRTS(N)  *  DSGRT ( N*1.00) 

5  DIRT(N)  =  l.DO/DRTS(N) 

10  CCTHET  =  DCCS I THET  A ) 

SITHET  =  DSIN(THETA) 

IF( IFLAG.NE. 1. AND. THETA. NE.O. DO) SITHI  =  I.C0/S1THCT 
COMPUTE  THE  LEGENDRE  FUNCTIONS  . 


15 


16 


RLNN ( 1 )  =  1 . DO 

RLNN ( 2 )  =  SITHET  *ORT  S (  3) 

DU  15  N l  =  3, Ml 
N  =  Nl-l 
N2  =  2*N 

RLMN(Nl)  =  DRTS 1 N2  +1 ) *D I RT  C N2 ) *S I  THE T*RLNN ( N l - 1 ) 
IFIM.GT.l)  GO  TO  20 
IF(M.EQ.O)  GO  TO  16 
RLEG (21  =  RLNN ( 2 ) 

RLEG!  3)  *  ORTS(5)*COTHET*RLEG(2) 

GO  TG  20 

RLEG ( 1 )  *  1.00 

RLEG  (  2)  *  CGTHET*'JRTS(  3) 
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20  CONTINUE 

RLEG (Ml  )  =  RLNN(MI) 

RLEG ( M2 )  *  DRTS(  Ml*2+l  ) *CUTHCT*RL EG( M 1 J 
DO  30  N I  =  M3 1 NMX l 

N  =  Nl-l 

IF(M.EQ.0.AN0.N.LT.2  .OR.  M.  ED  .  1  . AND.N.  L  T  .  3  )  GU  TO  30 
N2  =  2  *N 

RLEGCNl  )  =  DRTSl  N2  +  1)  *01  RMN  +  M  )*DIRT(N-M  )  *{DRTS  (M2-1 )  *CCTHET* 

2  i\LEG(  N  I— I )— DRTS  (N+M—  l J  *DRTS ( N-M- L l *0 l RT  l  N2-3 )  *RLEG(Nl-2)  ) 

GO  TO  30 
30  CONTINUE 

IF( IFLAG.Eu. 1)  RETURN 
IFtSITHET .tQ.O.DO)  wRITE(6,99I 

99  FORMAT!//*  ***  LEGFON  DOES  NOT  COMPUTE  DERIVATIVES  AT  THE  POLES 
2  *****.******$*****•  // ) 

IFl S ITHET .EG •0*00)  RETURN 


COMPUTE  ALL  THE  DERIVATIVES  OF  THE  LEGENDRE  FUNCTIONS. 


40 


2 


60 


RLNN(i)  =0.00 
KLN  =  RLNNI2J 
RLNNI2)  =  DRTSl 3 
DO  40  Ni  = 

N  =  Nl-l 
N2  =  2*N 
KLN 1  =  RLNN(Nl) 

R  L  NN (Nil  =  DRTS ( 
RLN  =  RLN1 
CONTINUE 


DLEG(MI)  = 

K  LNN  ( 

DO  60 

M  = 

N  =  Nl-l 

N2  =  N*2 

DLEG(Nl)  = 

S  I  THI 

DKT  S ( N2> l ) 

♦  DIRT 

CONTINUE 

RETURN 

END 

)*CUTHET 
3,  Ml 


N2*l)*0lRT(N2  )*(SITHET*RLNMN)+CDTHFT*RLNJ 


Ml  ) 

M2  ,  NMX  l 


*(  N  ♦RLEGINl  I*COTHET-DRrS(N-M)*ORTS(.\+M 
(  N2-l)*RLEGtN)  ) 


)* 


SUBROUTINE  kAOGRA( FL AT , FLON, H T ,KL AT , GR ,R E ) 

IMPLICIT  REAL*d( A-H,C-Z) 

C  THIS  SUBROUTINE  COMPUTES  GEOCEMTkIC  DISTANCE  TO  THE  POINT, 

C  THE  GEOCENTRIC  LATITUDE, AND 

C  AN  APPROXIMATE  VALUE  OF  NORMAL  GRAVITY  AT  THE  POINT  BASED 
C  CN  CONSTANTS  CF  THE  GEODETIC  REFERENCE  SYSTEM  198G 

DATA  AL/o 373 137 . 00/ , E2 /. 00669433002 2900/ ,R AD/ 57 .293  7  79  61 i0  32  32D0/ 

REAL*3  N 

FLATR=FLAT/R AD 

FLONR  =FLON/ R AD 

T1=0SIN(FLATR)**2 

N=AE/DS£R  T  (  l  .-E2  *T  I  J 

T2  =  (N«-HT)  *DCOS(  FLATR) 

X=T2*DC03< FLONP J 
Y  =  T2*0SIN(FLCNR  ) 

Z=l N*l l .-E2) +HT  >  *OSIN( FLATR) 

N  =  AE/DSQRT( 1 ,-62*Tl ) 

C  COMPUTE  THE  GEOCENTRIC  RADIUS 
RE=OScRn  X**2+Y**2  +  Z**2l 
C  COMPUTE  THE  GEOCENTRIC  LATITUDE 

RLAT=0ATAf,{  Z/DSORT  (  X**2+Y**2)  ) 

C  COMPUTE  NORMAL  GRA V  I T Y : UN I TS  ARE  M/SEC**2 

GR=9. 730 326  7  7 1500*  (  l  •  *-.00 19 31851  3 53DG*Tl )  /O SQKT  (  l.-E2*TI  I 
C  CORRECT  FOR  ELEVATION  IN  AN  APPROXIMATE  m AY 
GR=GR— HT*0 .30360—5 
RETURN 
END 
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